// Judicial Decisions, Police Officer Uncertainty, and the Escalation of Force //
// Courtenay R. Monroe, Sophia Hatz and Kristine Eck//
// MAIN ANALYSIS//
// December 2023 //
// note: this code draws on replication for de Chaisemartin and Hoetfouille 2020


ssc install did_multiplegt, replace
ssc install estout, replace

clear

use "JD_replication_data.dta"



// set panel data set

xtset state year

// check that making diff variable in Stata is the same as the diff variable from R
g restrictive_aele_diff_st = d.restrictive_aele 
tab restrictive_aele_diff_st 
tab restrictive_aele_diff_st  year
* from R:
tab restrictive_aele_diff 
tab restrictive_aele_diff year

// make dummy variables for year
tab year, gen(d)
// make dummy variables for state
tab state, gen(st)



* table 2
qui{

// m1
set seed 1
did_multiplegt police_shoot_fe state year restrictive_aele, breps(200) cluster(state)
graph drop _all 
scalar did_m1=e(effect_0)
scalar se_did_m1 = e(se_effect_0)
scalar N_did_m1 = e(N_effect_0)
scalar t_stat_m1 = e(effect_0)/e(se_effect_0)
scalar p_val_m1 = 2*normal(-abs(t_stat_m1))


// m1b
set seed 1
did_multiplegt police_shoot_fe state year restrictive_aele, controls(homicide_rate_ucr property_crime_rate_ucr) breps(200) cluster(state) 
graph drop _all 
scalar did_m1b=e(effect_0)
scalar se_did_m1b = e(se_effect_0)
scalar N_did_m1b = e(N_effect_0)
scalar t_stat_m1b = e(effect_0)/e(se_effect_0)
scalar p_val_m1b = 2*normal(-abs(t_stat_m1b))

// m2
set seed 1
did_multiplegt police_shoot_fe state year restrictive_count_aele, recat_treatment(restrictive_count_aele_recat) breps(200) cluster(state)
graph drop _all 
scalar did_m2=e(effect_0)
scalar se_did_m2 = e(se_effect_0)
scalar N_did_m2 = e(N_effect_0)
scalar t_stat_m2 = e(effect_0)/e(se_effect_0)
scalar p_val_m2 = 2*normal(-abs(t_stat_m2))

// m2b
set seed 1
did_multiplegt police_shoot_fe state year restrictive_count_aele, recat_treatment(restrictive_count_aele_recat) controls(homicide_rate_ucr property_crime_rate_ucr) breps(200) cluster(state)
graph drop _all 
scalar did_m2b=e(effect_0)
scalar se_did_m2b = e(se_effect_0)
scalar N_did_m2b = e(N_effect_0)
scalar t_stat_m2b = e(effect_0)/e(se_effect_0)
scalar p_val_m2b = 2*normal(-abs(t_stat_m2b))



matrix res1 = (did_m1, se_did_m1,t_stat_m1,p_val_m1, N_did_m1\ ///
			  did_m1b, se_did_m1b,t_stat_m1b,p_val_m1b, N_did_m1b\ ///
			  did_m2, se_did_m2,t_stat_m2,p_val_m2, N_did_m2\ ///
			  did_m2b, se_did_m2b,t_stat_m2b,p_val_m2b, N_did_m2b)
			  
matrix rownames res1=m1 m1b m2 m2b 
matrix colnames res1=estimate se tstat pval N
}

// 
matrix list res1

putexcel set "table2", replace
putexcel A1=matrix(res1), names


* pre-treatment trends in police violence
//m1
set seed 1
did_multiplegt police_shoot_fe state year restrictive_aele, placebo(5) breps(200) cluster(state)
scalar did_m1=e(effect_0)
scalar se_did_m1 = e(se_effect_0)
scalar pl1_m1=e(placebo_1)
scalar se_pl1_m1 = e(se_placebo_1)
scalar pl2_m1=e(placebo_2)
scalar se_pl2_m1 = e(se_placebo_2)
scalar pl3_m1=e(placebo_3)
scalar se_pl3_m1 = e(se_placebo_3)
scalar pl4_m1=e(placebo_4)
scalar se_pl4_m1 = e(se_placebo_4)
scalar pl5_m1=e(placebo_5)
scalar se_pl5_m1 = e(se_placebo_5)

matrix m1p = (did_m1, se_did_m1\pl1_m1,se_pl1_m1\ ///
pl2_m1,se_pl2_m1\pl3_m1,se_pl3_m1\ ///
pl4_m1,se_pl4_m1\pl5_m1,se_pl5_m1)
			  
matrix rownames m1p=t0 t1 t2 t3 t4 t5 
matrix colnames m1p=estimate se 
matrix list m1p

putexcel set "m1placebo", replace
putexcel A1=matrix(m1p), names


// m2
set seed 1
did_multiplegt police_shoot_fe state year restrictive_count_aele, recat_treatment(restrictive_count_aele_recat) placebo(5) breps(200) cluster(state)
scalar did_m2=e(effect_0)
scalar se_did_m2 = e(se_effect_0)
scalar pl1_m2=e(placebo_1)
scalar se_pl1_m2 = e(se_placebo_1)
scalar pl2_m2=e(placebo_2)
scalar se_pl2_m2 = e(se_placebo_2)
scalar pl3_m2=e(placebo_3)
scalar se_pl3_m2 = e(se_placebo_3)
scalar pl4_m2=e(placebo_4)
scalar se_pl4_m2 = e(se_placebo_4)
scalar pl5_m2=e(placebo_5)
scalar se_pl5_m2 = e(se_placebo_5)

matrix m2p = (did_m2, se_did_m2\pl1_m2,se_pl1_m2\ ///
pl2_m2,se_pl2_m2\pl3_m2,se_pl3_m2\ ///
pl4_m2,se_pl4_m2\pl5_m2,se_pl5_m2)
			  
matrix rownames m2p=t0 t1 t2 t3 t4 t5 
matrix colnames m2p=estimate se 
matrix list m2p

putexcel set "m2placebo", replace
putexcel A1=matrix(m2p), names







// FD reg
* m1
reg police_shoot_fe_diff restrictive_aele_diff d2-d19, cluster(state)
scalar b_fd_m1=_b[restrictive_aele_diff]
*m1b
reg police_shoot_fe_diff restrictive_aele_diff homicide_rate_ucr_diff property_crime_rate_ucr_diff  unemp_rate_ukcpr_diff poverty_rate_ukcpr_diff population_ukcpr_log_diff d2-d19, cluster(state)

* m2
reg police_shoot_fe_diff restrictive_count_aele_diff d2-d19, cluster(state)
scalar b_fd_m2=_b[restrictive_count_aele_diff]
* m2b
reg police_shoot_fe_diff restrictive_count_aele_diff homicide_rate_ucr_diff property_crime_rate_ucr_diff  unemp_rate_ukcpr_diff poverty_rate_ukcpr_diff population_ukcpr_log_diff d2-d19, cluster(state)


// TWFE reg
* m1
reg police_shoot_fe restrictive_aele d2-d19 st2-st51, cluster(state) 
scalar b_fe_m1=_b[restrictive_aele]
*m1b
reg police_shoot_fe restrictive_aele homicide_rate_ucr property_crime_rate_ucr  unemp_rate_ukcpr  poverty_rate_ukcpr  population_ukcpr_log  d2-d19 st2-st51, cluster(state) 

* m2
reg police_shoot_fe restrictive_count_aele d2-d19 st2-st51, cluster(state) 
scalar b_fe_m2=_b[restrictive_count_aele]
* m2b
reg police_shoot_fe restrictive_count_aele homicide_rate_ucr property_crime_rate_ucr  unemp_rate_ukcpr  poverty_rate_ukcpr  population_ukcpr_log  d2-d19 st2-st51, cluster(state) 


// Weights attached to the FE and FD regressions 


* FD weights
twowayfeweights police_shoot_fe_diff state year restrictive_aele_diff restrictive_aele, type(fdTR) 
display 48/(48+191)
twowayfeweights police_shoot_fe_diff state year restrictive_count_aele_diff restrictive_aele, type(fdTR) 
display 68/(68+171)

* correlation b/w weights and year
twowayfeweights police_shoot_fe_diff state year restrictive_aele_diff restrictive_aele, type(fdTR) test_random_weights(year) 
twowayfeweights police_shoot_fe_diff state year restrictive_aele_diff restrictive_count_aele, type(fdTR) test_random_weights(year) 


* FE weights 
twowayfeweights police_shoot_fe state year restrictive_aele, type(feTR)

twowayfeweights police_shoot_fe state year restrictive_count_aele, type(feTR)
display 43/(43+196)

* correlation b/w weights and year
twowayfeweights police_shoot_fe state year restrictive_aele, type(feTR) test_random_weights(year) 
twowayfeweights police_shoot_fe state year restrictive_count_aele, type(feTR) test_random_weights(year) 




// Tests that beta_FE, beta_FD and DID_M are equal (see 3rd paragraph p.23)


scalar diff1=b_fe_m1-b_fd_m1
scalar diff2=b_fe_m1-did_m1	
scalar diff3=b_fd_m1-did_m1
matrix B=(0,0,0)

*m1
set seed 1
qui{
forvalue i=1/200{
preserve
bsample, cluster(state)
xtreg police_shoot_fe d2-d19 restrictive_aele, fe cluster(state)
scalar b_fe_m1=_b[restrictive_aele]
reg police_shoot_fe_diff restrictive_aele_diff d2-d19, cluster(state)
scalar b_fd_m1=_b[restrictive_aele_diff]
did_multiplegt police_shoot_fe state year restrictive_aele, cluster(state)
scalar did_m1=e(effect_0)
matrix B=B\(b_fe_m1-b_fd_m1,b_fe_m1-did_m1,b_fd_m1-did_m1)
restore
}
}
preserve
drop _all
svmat B
drop if _n==1
forvalue i=1/3{
qui{
su B`i'
}
if `i' == 1 { 
	di "The t-stat for the difference b/w TWFE and FD is "
}
else if `i' == 2 {
	di "The t-stat for the difference b/w TWFE and DID_M is "
}
else {
	di "The t-stat for the difference b/w FD and DID_M is "
}
di diff`i'/r(sd)
}
restore


*m2
set seed 1
qui{
forvalue i=1/200{
preserve
bsample, cluster(state)
xtreg police_shoot_fe d2-d19 restrictive_count_aele_recat, fe cluster(state)
scalar b_fe_m1=_b[restrictive_count_aele_recat]
reg police_shoot_fe_diff restrictive_count_aele_diff d2-d19, cluster(state)
scalar b_fd_m1=_b[restrictive_count_aele_diff]
did_multiplegt police_shoot_fe state year restrictive_count_aele, recat_treatment(restrictive_count_aele_recat) cluster(state)
scalar did_m1=e(effect_0)
matrix B=B\(b_fe_m1-b_fd_m1,b_fe_m1-did_m1,b_fd_m1-did_m1)
restore
}
}


preserve
drop _all
svmat B
drop if _n==1
forvalue i=1/3{
qui{
su B`i'
}
if `i' == 1 { 
	di "The t-stat for the difference b/w TWFE and FD is "
}
else if `i' == 2 {
	di "The t-stat for the difference b/w TWFE and DID_M is "
}
else {
	di "The t-stat for the difference b/w FD and DID_M is "
}
di diff`i'/r(sd)
}

restore

* alternative operationalisation (reported in appendix)

* table 4
qui{

// m3
set seed 1
did_multiplegt police_shoot_fe state year status_quo, breps(200) cluster(state)
graph drop _all 
scalar did_m3=e(effect_0)
scalar se_did_m3 = e(se_effect_0)
scalar N_did_m3 = e(N_effect_0)
scalar t_stat_m3 = e(effect_0)/e(se_effect_0)
scalar p_val_m3 = 2*normal(-abs(t_stat_m3))


// m3b
set seed 1
did_multiplegt police_shoot_fe state year status_quo, controls(homicide_rate_ucr property_crime_rate_ucr) breps(200) cluster(state)
graph drop _all 
scalar did_m3b=e(effect_0)
scalar se_did_m3b = e(se_effect_0)
scalar N_did_m3b = e(N_effect_0)
scalar t_stat_m3b = e(effect_0)/e(se_effect_0)
scalar p_val_m3b = 2*normal(-abs(t_stat_m3b))


matrix res2 = (did_m3, se_did_m3,t_stat_m3,p_val_m3, N_did_m3\ ///
			  did_m3b, se_did_m3b,t_stat_m3b,p_val_m3b, N_did_m3b)
			  
matrix rownames res2=m3 m3b 
matrix colnames res2=estimate se tstat pval N
}

// 
matrix list res2
putexcel set "table4", replace
putexcel A1=matrix(res2), names

